Towards an efFicient multiscale modeling of low-dimensional reactive systems: 
study of numerical closure procedures 

Giacomo Mazzi/'[^ Yannick De Decker,^ and Giovanni Samaey-"^ 

Department of Computer Science, KU Leuven, Belgium 
^^Interdisciplinary Center for Nonlinear Phenomena and Complex Systems (CENOLI) , Universite Libre de Bruxelles, 
Belgium 

In this paper, we present a study on how to develop an efficient multiscale simulation strategy for the 
dynamics of chemically active systems on low-dimensional supports. Such reactions are encountered in a 
wide variety of situations, ranging from heterogeneous catalysis to electrochemical or (membrane) biological 
processes, to cite a few. We analyzed in this context different techniques within the framework of an important 
multiscale approach known as the equation free method (EFM) , which "bridges the multiscale gap" by building 
microscopic configurations using macroscopic-level information only. We hereby considered two simple reactive 
processes on a one-dimensional lattice, the simplicity of which allowed for an in-depth understanding of the 
parameters controlling the efficiency of this approach. We demonstrate in particular that it is not enough to 
base the EFM on the time evolution of the average concentrations of particles on the lattice, but that the time 
evolution of clusters of particles has to be included as well. We also show how important it is for the accuracy 
of this method to carefully choose the procedure with which microscopic states are constructed, starting 
from the measured macroscopic quantities. As we also demonstrate that some errors cannot be corrected by 
increasing the number of observed macroscopic variables, this work points towards which procedures should 
be used in order to generate efficient and reliable multiscale simulations of these systems. 



INTRODUCTION 

Chemically reacting systems exhibit behaviors that 
may be described at different levels of modeling. At 
a fine microscopic level, the system consists of a large 
number of atoms and molecules (particles) of different 
species, and reactions at this level can be seen as discrete 
events that occur with some probability when two or 
more particles are close enough to interact. At a coarser, 
macroscopic level one only considers the concentration of 
each species, whose evolution is modeled via a system of 
ordinary or partial differential equations (depending on 
whether spatial inhoniogeneities are taken into account). 
Whereas the microscopic level of description might allow 
modeling inter-particle interactions in great detail, more 
macroscopic models are clearly desirable for at least two 
reasons. First, fine-scale simulation of systems consisting 
of many particles may be prohibitively expensive, espe- 
cially when the chemical reactions span a wide range of 
time-scales. Second, even if feasible, one is usually only 
interested in the evolution of some coarse-level, macro- 
scopic quantities, such as species concentrations. 

When considering only macroscopic quantities, the dy- 
namics of the system is often described using a mean field 
approximation (MFA), which expresses the law of mass 
action for the species concentrations; the und erlyin g as- 
sumption being that the system is well stirrecP^^. Un- 
fortunately, in many systems of practical interest, this re- 
quirement is not fulfilled. A typical example is the case of 
catalytic reactions on surfaces, in which the particles are 
located on a spatial structure with low dimensionality. 
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such as a lattice. Then, the location of the potentially re- 
acting particles becomes important, as chemical reactions 
can only fire if the reactants are physically close to each 
other and mixing is so inefficient that local fluctuations 
of composition are not easily wiped out. A possible ap- 
proach is to expand the set of macroscopic state variables, 
taking into account particle correlations. As for the mean 
field equations, closed equations for the pair correlations 
can only be obtained using additional assumptions. Sev- 
eral analytical closure approximations have been devel- 
oped such as the quasi-chemical approximatiorP^, addi- 
tional references are given in Section [II A| A drawback of 
analytic closure approximations is that they may depend 
strongly on the system considered and when, besides pair 
correlations, other additional variables are needed, the 
equations quickly become intractable. 

When no closed macroscopic model with sufficient ac- 
curacy is available, one needs to resort to an individual- 
based stochastic particle simulation. In this case, the 
microscopic state of the system is updated via a kinetic 
Monte Carlo (kMC) algorithm according to a given set 
of reaction laws, such as the stochastic simulation algo- 
rithm (SSA)P, with adaptations to also take into account 
the underlying lattice structure in the aforementioned 
case of surface reactions. These simulations are often 
time consuming and several methods have been devel- 
oped that take advantage of disparate reaction rates to 
accelerate them. We mention, for instance, tau-leaping'^, 
R-leapin^, and nested SSA'^. However, most of these 
methods are not easily implemented in lattice kMC, un- 
less other approximations are also taken into account; see 
[25j for an example which makes use of a type of coarse- 
graining. For a general review of kinetic Monte Carlo 
methods in this context, see [3]. 

There exists, however, an alternative way of extract- 
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ing the macroscopic dynamics. Computational strategies 
have recently been developed that are capable of acceler- 
ating individual-based particle simulations by introduc- 
ing on-the-fly numerical approximate macroscopic clo- 
sures: we mention equation-free^"^ and heterogeneous 
multiscale methods (HMM)^'^. In both approaches, the 
crucial steps are (i) to decide on the set of macroscopic 
state variables that are sufficient to determine future 
evolution of the system, and (ii) to define an operator 
that generates a microscopic state corresponding to a 
given macroscopic state. This second step is actually 
equivalent to prescribing the closure approximation; in 
the equation-free, resp. HMM, frameworks, this step is 
called lifting, resp. reconstruction. The converse opera- 
tor (mapping a microscopic state to a macroscopic one) 
is called restriction in both frameworks. Once these op- 
erators are available, one can construct a coarse time- 
stepper, which evolves the chosen macroscopic state vari- 
ables over some time interval by means of a three-step 
procedure: (i) lifting; (ii) microscopic simulation; and 
(iii) restriction. 

Equation-free methods have already been used to per- 
form n umeri cal bifurcation analyses of lattice chemical 
systempE^l. in [T7], one studies a CO oxidation process 
on a lattice for which a bistability observed in simula- 
tions is not fully reproduced by a mean-field or a quasi- 
chemical approximation. Using a coarse time-stepper in 
which the macroscopic state is given by the concentra- 
tions of all species, however, the observed macroscopic 
bifurcation results can be completely reproduced and an- 
alyzed. This result is the consequence of two fundamental 
underlying properties of the system: (i) the coverages by 
themselves are indeed sufficient to predict macroscopic 
evolution of the system; (ii) given the global concentra- 
tions, one can construct a lifting operator that generates 
appropriate initial conditions for individual lattice sites. 
In [15] NO reduction by CO over a Pt(lOO) surface is 
also analyzed with the same approach. In both cases, 
the lifting equilibrates the lattice by exploiting a separa- 
tion in time scales between (fast) diffusion on the lattice 
and (slower) chemical reactions. As a result, all high or- 
der spatial correlations are slaved to the concentrations. 
This separation of time scales allows for what the authors 
call a property of "healing" : incorrectly initialized higher 
order correlations quickly relax to a correct functional of 
the coverages. 

The present paper assesses the efficiency of lift- 
ing/reconstruction techniques in the case of chemical re- 
actions on lattices; the procedure we propose, however, 
could be appHed to many multiscale models, see e.g. [55] 
for a similar study in the context of polymeric fluids. In 
particular, we investigate systems that are subject to a 
set of chemical reactions in a crowded environment that 
makes diffusion ineffective. Such situation puts in danger 
the aforementioned diffusion-induced separation of time 
scales. 

As already pointed out in [T7], two problems are ex- 
pected to arise in this context. First, we will more 



than probably need to augment the set of macroscopic 
state variables, as an approximate closure on the level 
of coverages alone will generally not suffice. Second, in 
the absence of a clear time-scale separation, artifacts in- 
troduced during the lifting step cannot be expected to 
"heal" completely. Additionally, we are not only inter- 
ested in capturing the correct steady state behavior in 
the present paper, but also in recovering the (transient) 
macroscopic dynamics. In this view, the main contribu- 
tion of the present paper is twofold: 



• From an algorithmic point of view, we develop a nu- 
merical closure strategy to recover the macroscopic 
(transient) dynamics of two specific chemical sys- 
tems with nonlinear reactions on a one-dimensional 
lattice, defined in Section |lj We introduce a hier- 
archy of macroscopic state variables for which the 
implementation of a consistent lifting operator is 
feasible (Section |ll]) and consider several options 
to initialize the remaining degrees of freedom (Sec- 
tion Iml. 



• From a computational physics point of view, we 
analyze the effect of the different lifting strategies 
(corresponding to different closure approximations) 
on the macroscopic evolution (Section IV). This 



study reveals that both the number and type of 
macroscopic state variables and the choice of the 
lifting operator (for given macroscopic state vari- 
ables) can have a significant effect on the observed 
macroscopic evolution and hence, on the accuracy 
of the multiscale approach. 



We analyzed this on the case of two chemical systems 
with nonlinear reactions on a one-dimensional lattice. 
Our choice of one-dimensional systems is justified not 
only by a corresponding simple analysis, but also because 
this setting has the additional advantage that effects due 
to the low dimensionality are more pronounced. Conse- 
quently, if a procedure provides accurate results in one 
space dimension, we expect it to also perform well in 
higher dimensions. 

The paper is organized as follows: in Section [T| we 
present the model problems considered in this work. We 
describe how these systems are modeled, as well as how 
they are simulated numerically. In Section [lT| we intro- 
duce the multiscale framework, including the hierarchy of 



macroscopic state variables (Section II A), the different 
lifting and restriction operators (Se ction |IIB ) , and the 
numerical closure algorithm (Section II C I. In Section III 
we describe the different lifting procedures. The numeri 
cal results we have obtained on two model problems are 
presented in SectiorjlVj We conclude in Section |Vj where 
we also outline open questions and directions for future 
research. 
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I. MODEL PROBLEMS AND COARSE-SCALE 
APPROXIMATIONS 

A. Chemical reactions on one-dimensional lattices 

We consider a one-dimensional spatial domain with 
N lattice sites and periodic boundary conditions (i.e., a 
ring). Each lattice site is occupied by exactly one parti- 
cle. For simplicity, we consider only two particle species 
X and A. The dynamics of the system is governed by 
a set of / reactions, with associated reaction rates ki 
(1 < i < /). Reactions can only occur if the reacting 
particles are in adjacent sites. We adopt a slight modi- 
fication of the n otations that were introduced previously 
in the literatur^SIIlIlil .^^g define the microscopic state 
at time t as a{t) = {<Ti{t), a-2{t), . . . , a-]y{t)} where each 
variable (T„(t) = ±1, depending on whether the lattice 
site n is occupied by A (cr„(t) = +1) or X (cr„(t) = —1). 
We also define by 



ait) 



(1) 



1 



the coverage of species A, and correspondingly x 
the coverage of species X. 

The system state evolves in time as a consequence of 
chemical reactions. Throughout this text, we will con- 
sider two particular example systems: 

Example I.l (Schlogl model). Consider the following 
two-species system, in which two reversible reactions are 
possible, 



XAX 
X 



XXX, 
A, 



(2) 



where the notation in the first line is chosen to empha- 
size that only elements surrounded by two X may react. 
This type of system is a two-species restriction of the 
popular Schlogl model (usually called Schlogl's second 
modefl^, and it is widely used as prototype for study of 
non-equilibrium systems on constra ined c onfigurations, 
such as catalytic reactions on surfaceJ^^EH. Similar equa- 
tions have been proposed to describe disease spreading in 
epidemic models'*. 

Example 1.2 (Simple trimolecular reaction). We may 
simplify the above prototypical example even further by 
only retaining the trimolecular reaction, i.e.. 



XAX 



XXX, 



(3) 



where we assume fci = /c2 = 1. This model 
contains the simplest mechanism to obtain a macro- 
scopic behavi or th at non-trivially depends on the local 
configuration421124] Note that this example arises when 
considering Schlogl's second model in the limit where 
^3 = ^4 = 0. 



B. Stochastic simulation 

A common way to simulate chemical reactions on a 
lattice is via stochastic simulation, in which one com- 
putes one realization of the stochastic process as a func- 
tion of time. As in most of the earlier literature for 
and Example O UJMi^ .^g apply a kinetic 
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Example 

Monte Carlo (kMC) method with an accep t- rej ect strat- 
egy (a so-called null-event kMC, see, e.g.^^^. As for 
the classical stochastic simulation algorithm (SSA|^ for 
well-stirred systems (which is rejection- free), some mod- 
ifications are required, as the reaction probabilities not 
only depend on the selected reaction but also on the cho- 
sen lattice site. 

Any kMC method depends on the definition of a (pos- 
sibly time-dependent) propensity function Vi(t) per re- 
action (1 < « < /); the quantity Ti[t)dt represents the 
probability of reaction i occurring in the time interval 
[t,t + dt\. The main modification to accommodate for 
the simulation of reactions on low-dimensional lattices is 
that the propensity functions, which are spatially inde- 
pendent for well-stirred systems, now also depend on the 
reaction site n = 1, . . . ,N . Therefore, at each moment in 
time, we have N x I propensity functions, which we de- 
note as Ti^n, with 1 < n < N. Moreover, the propensity 
functions do not depend on the total coverages of each 
of the species, but on t he states of the individual lattice 

may be derived from ^ as 

l + a„it) (l-cr„_i(f)) 



sites. For Example 
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XAX 



XXX 



k2. 



X 



A 



XXX, 


Ti 


n(i) 


XAX, 




n(i) 


^A, 




n(i) 


^X, 


r4 


n(i) 



2 

2 

-cr„(t). 



lit)) (1 - a„ 



(4) 
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We further introduce the total propensity for a lattice 
site n as the sum of all reaction propensities for that 
lattice site, 



I 



where Ttc,t,nit)dt represents the probability that any re- 
action will occur at the site n in the time interval [t, t+dt]. 
The maximal total propensity over all lattice sites is de- 
fined as rmax(i) = inaxn rtot,n(i). Finally, we will also 
use the total propensity of a reaction i, which is defined 
as the sum of the reaction propensities for that reaction 
over all lattice sites 

JV 

r^,tot(t) = Er,^„(t). 

n=l 

One time step of the algorithm then proceeds as follow^ 

1. A lattice site m is randomly chosen (from a dis- 
crete uniform distribution on {1, . . . , N}), and all 
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propensity functions Ti^mit) for the reactions are 
calculated at that site, and renormalized to obtain 
fi,m(t) = ri^,„(t)/rniax(^)- Notc that, as a conse- 
quence of the renormalization, we have the prop- 
erty 



< 1. 



(5) 



2. Next, a reaction is randomly chosen according to 
the propensity functions. First generate a random 
number ^ from the uniform distribution on [0,1]. 
The chosen reaction is then given as j such that 



(6) 



i=l 



i=l 



When rtot,m(i) < 1, we might have that 
rtot,m(i)/rmax(i) < C < 1- In that case, there is 
no reaction j < I for which the above equation is 
satisfied; this situation corresponds to a null-event, 
and no reaction is performed. By convention, we 
set j = / + 1 in this case. 



3. Compute the time increment 



i/r„tot(0, 





1 < J < / 

J = 1 + 1 



and advance the system clock, t = t + Stj{t). 

The last step differs from SSA, as there the size of Stj is 
evaluated using another random generated number. 



Remark I.l. In the case of Example |I.2[ this algorithm 
may be simplified considerably. We can then pick ran- 
domly a lattice site m (again according to a discrete uni- 
form distribution). If cr„i-i(i) = (Tm+i(t) = —1, then we 
change the state of that site a„Xi + <^i) — For 
this example, bt is constant as a function of time, and 
independent of the chosen reaction. 



C. Master equation and approximation of macroscopic 
dynamics 

Instead of performing a stochastic simulation that gov- 
erns the evolution the configuration for an individual re- 
alization of the system, one can also describe the deter- 
ministic evolution of the probability distribution of the re- 
alizations over the configuration space, which is governed 
by the master equation. While the computational cost of 
this approach is prohibitive for practical purposes, this 
description is useful to derive approximate macroscopic 
models. 

Let us introduce the probability distribution 
7'(cr, tjSE^, and corresponding master equation 



N 



Wn{(T)V{(T,t)], (7) 



in which a = {ai, . . . , cr„, . . . , ctat}, and each o'" = 
{fJi, . . . , (T„_i, — (T„, (T„+i, . . . , ctat} differs from cr in ex- 
actly one lattice site. The rates u'„(cr) are the probabil- 
ity per unit time that the n-th particle fiips from tT„ to 
— (T„, given that the current state of the whole system 
is (T. For the systems we consider in the present paper, 
the form of Wn depends on the different reaction laws 
and on the nearest neighborhood {(t„_i, ovi+i}. For in- 
stance, given the reaction laws of Example |I.H we get the 
reaction probabilities per sit^^o)^ 



w„ (a) 



+ 



ki + k2 ki 
2 ^ 
ki + ka 



2 

k4 - ks 



1 — 1 — CTn + l 



(8) 



The Wn are also related with the propensity functions Q 
as w„(cr(t)) = Ttot,n{t)- Note that the first two terms, 
which refer to the first reaction, depend on the state of 
the nearest neighbors, in contrast to the last two terms, 
which refer to the second reaction. 

For a majority of systems, the master equation ([7| is 
too complicated to be solved directly. However, when 
introducing suitable approximations, it is possible to ex- 
tract from it equations of motions for appropriately cho- 
sen macroscopic quantities (such as coverages). The 
simplest approximate macroscopic model is the so-called 
mean field approximation, where, by substituting all the 
high order correlations with an average "mean field" in- 
teraction, one writes a set of ordinary differential equa- 
tions for the coverage j^^ * ^° l that can be simulated with 
drastically reduced computational complexity. 

Assuming translational invariance, it was shown iipSlini 
that the mean-field evolution equation for q{t) = {an{t)) 
is here given by 

^ = -^(fci - fcs) + (^3 - ki) + f|(fci - Sfcs) + k3 + ki]q 



(9) 



which may be written in terms of one of the coverages 
(here x = (1 — q)/2), and reads 



dx 
'dt 



-(fci -I- k2)x^ + kix^ - (fcg + ki)x + ki. (10) 



If we set k^ = ki = 0, we get the equations for Example 

in 



dx 
'dt 



= — fcia;'^ -|- k2ax'^ 



(11) 



The mean-field approximation is based on neglecting 
of all pair (and higher order) correlations. It is therefore 
expected to perform poorly in the systems of interest in 
this paper: in fact, this level of description even yields 
equilibrium values that differ substantially from those ob- 
tained by the microscopic dynamics Q^^l neverthe- 
less represents a suitable starting point for more involved 
approximation strategies. For instance, it is possible to 
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also consider pair correlations via the Ursell expansiorP. 
This approach consists in a development of pair, triplet, 
quadruplet, etc. correlation functions from which one 
can systematically extract finer and finer descriptions as 
higher order terms are kept. Another approach consists 
in expressing high order moments directly as functional 
of the lower ones. This is the case for the popular Kirk- 
wood approximation scheme^ or the (less well-known) 
Bethe-type ansatJ^^l xhe efficiency of several of these 
approaches for the models under consideration here are 
discussed in [12] where it was shown that 3-points corre- 
lations play as an important role as pair correlations in 
the dynamics. In any case, despite the simplicity of the 
considered model, the derivation and analysis of higher 
order approximations already revealed to be quite tedious 
and intricate. This justifies the development of numerical 
closure techniques, which we will discuss now. 



II. NUMERICAL CLOSURE ALGORITHM 

We now describe the numerical closure algorithm that 
is the focus of the present paper. We first define a hierar- 
chy of macroscopic state variables (Section II A I. Subse- 



quently, we introduce numerical lifting and restriction 
procedures to map macroscopic state variables to mi- 
croscopic realizations, and vice versa (Section II B I. Fi- 



nally, we propose a coarse time-stepper that mimics the 
evolu tion o f the corresponding unavailable macroscopic 
modePES. 



A. A hierarchy of macroscopic state variables 

As in the mean-field approximation, we first introduce 
a macroscopic state variable a{t), see equation Q, that 
represents the coverage of species A; recall that the cov- 
erage of species X is obtained as x = 1 — a. In addi- 
tion to the coverage, we introduce a hierarchy of macro- 
scopic state variables. Each macroscopic state variable 
1 < I < La, (respectively Ml^, I < I < Lx) is 
defined as the number of clusters of particles A (respec- 
tively X) of exactly length Z up to a maximum length of 
La (respectively Lx)- Explicitly, these additional state 
variables can be written as 



^ N l-l 



n=l 
N 



k=0 

-1 



^ ^{cTn-l + 1)(ct„+( + 1) J|(cr„+fc - 1), 



n=l 



k=0 



(12) 

where periodic boundary condition on the index n are 
assumed. The complete set of macroscopic state variables 
is then denoted as 

= {a, MA , . . . , Mi , . . . , }. 



Remark II. 1. We stress that the above definition is 
slightly different from other definitions of macroscopic 
state variables based on cluster counting already pre- 
sented in the literature, the most notable being the Kirk- 
wood approximation^. For instance, the variable 
only takes into account single A's (surrounded on both 
sides by X). Also, a cluster of length L^a,x} = 3 is not 
counted simultaneously as a two clusters of length 2, as 
is done in the Kirkwood approximation; we only count 
structures of the form AXXA to determine the number 
of clusters of length 2. To avoid potential confusion, we 
introduce the notation 



{XXX)^ = 



1 



2'H 



N l~l 



1), (13) 



to introduce the number of lattice sites occupied by 
species X, and also surrounded by species X, as it would 
appear in the Kirkwood approximation. 



As will become clear in Section II B this hierarchy of 



macroscopic state variables greatly simplifies the imple- 
mentation of a consistent lifting strategy, without intro- 
ducing additional complexity into the restriction. 



B. Lifting and restriction 

We introduce two operators that make the transition 
between microscopic and macroscopic state variables. 
We define a lifting operator, 



C : Ml, ^ <T, 



(14) 



which maps a macroscopic state to a realization of a chain 
of length N, and an associated restriction operator. 



(15) 



which maps a microscopic realization to the correspond- 
ing macroscopic state. The restriction operator is readily 
defined using formulas (12). The lifting step, however, is 
more involved. Besides the obvious condition that the 
lifting needs to be consistent, i.e., 

7^ o £ = Id, 

one also needs to ensure that microscopic configurations 
that result from lifting do not contain undesirable arti- 
facts, and actually represent realistic microscopic config- 
urations. For such purpose we introduce the concept of 
dynamical equivalence to compare different microscopic 
configurations. 

For the reaction laws ([2]), we say that two microscopic 
configurations cri and <T2 are equivalent if they have the 
same number of possibly reacting lattice sites for each of 
the possible reactions, i.e., if 



(XXX)^, = {XXX)^„ 

M\{rT,) = M\(cT2), 
a(cri) = a(cri), 



(16) 
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where the notation in the first equahty needs to be un- 



derstood in the Kirkwood sense, see equation (13). The 



last two conditions are satisfied for microscopic config- 
urations that have the same macroscopic state, for any 
Lx and La > 1, while the first condition is satisfied ex- 
actly only by conserving all the clusters (for instance, 
by choosing Lx = La = N). In this case, lifting cor- 
responds to a random shuffle of clusters: the number of 
cluster is conserved, but not their exact locations. So, 
using the definition of equivalence introduced above, one 
has 



*t = £(Ma,,jv(t)). 



(17) 



When decreasing the number of macroscopic state vari- 
ables (by choosing smaller Lx and La)^ the numbers of 
longer clusters will not be preserved during lifting. In- 
stead, we only conserve the total number of elements of 
species X and species A that are not contained within the 
counted clusters. These remainders can be computed as 



rA = aN -Y,IM\, 

1=1 

rx ={N - aN)-J2^M 



(18) 



I 

X- 



1=1 



During a lifting procedure, one will typically choose a 
limited amount of macroscopic variables from the avail- 
able set to construct the microstates. This reconstruc- 
tion can be done in many different ways, depending on 
how one treats the above-defined remainders to create a 
lattice state that is compatible with the macroscopic con- 
straints. We here emphasize that each such lifting policy 
actually corresponds to an implicit choice of closure. As 
it can clearly be seen now, lifting indeed corresponds to 
reconstructing higher order moments, starting from low 
order ones. The investigation of different choices for the 
lifting is exactly the focus of the present paper. In par- 
ticular, the main goal is to investigate which lifting op- 
erator recovers the macroscopic dynamics of the system 
with a minimal number of macroscopic variables. In Sec- 
tion |lll[ a number of choices will be presented. Numerical 
results using the resulting lifting operators are given in 
Section HVl 



for the microscopic model, consistent with the 
macroscopic state M^^^^^ at time t* . 

(ii) Simulation of the evolution of the microscopic state 
over a time interval [t* , t* + 5t] . 

(iii) Restriction, i. e. the observation (estimation) of the 
macroscopic state at t* + St: 

Mla,lx {t* + St) = n{(T{t* + St)). 

The coarse time-stepper can be used to accelerate 
simulation in time using coarse projective integration 
algorithmic^, or to perform a numerical bifurcation 
analysi^l. In this paper, we will use the coarse time- 
stepper to study how different choices for the numbers 
La and Lx of macroscopic state variable M^^ and 
of the specifics of the lifting (determining the distribution 
of the remainders and rx ) introduce a closure approx- 
imation (i.e., affect the macroscopic dynamics). To this 
end, we compare the time series of the macroscopic state 
variables M^^^^^ obtained via the numerical closure al- 
gorithm with those obtained by a complete microscopic 
simulation, Mi^_L^(i) = TZ{a{t)). A similar study has 
been performed for a model for polymeric fluidi^H 



III. LIFTING PROCEDURES 

Let us now detail the specific lifting operators that we 
propose in this paper. Each lifting operator consists of 
two steps: (i) positioning of the conserved clusters of 
species X and A, as specified by the macroscopic state 
variables M^^^/,^; and (ii) distribution of the remaining 
rA (respectively rx) elements of species A (respectively 
X). 

The first step is common to all developed strategies, 
limiting differences in lifting operators to the distribu- 
tion of the remainders. Each macroscopic state variable 
M\ (1 < Z < La), respectively Mj^ (1 < Z < Lx), 
refers to the number of clusters of length I for species 
A, respectively X. We gather these clusters in two sets: 
Ba = {A,...,A,AA,...AA,AAA,...,AAA,...}, and 
similarly Bx , and then select elements from Ba and Bx 
alternatingly until one of the two sets is empty. We fur- 
ther introduce the notation 



C. The numerical closure algorithm 

We are now ready to formulate the numerical closure 
algorithm that will be used to simulate time evolution of 
the macroscopic state variables. Given an initial condi- 
tion for the macroscopic state variables M(t*) at time i*, 
we define coarse time-stepper as the following three-step 
procedure: 

(i) Lifting, i.e. the creation of initial conditions 

cT{e) = c{mLAxAn), 



ruA 



\Ba\ 



La 
1=1 



mx 



\B 



X 



Lx 

1=1 



(19) 

for the cardinality of both sets. Since usually mA rnx, 
some elements from rA or rx will need to be used to 
ensure that we can place all clusters that are accounted 
for in the macroscopic state variables (see below) . In the 
remainder, we assume that mx > mA, the argument for 
the alternative case is identical. 

Given the definition of the macroscopic state variables, 
we need to group elements contained in the remainders 
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TA (respectively rx) in blocks longer than La (respec- 
tively Lx)- The lifting procedures described below differ 
exactly in the way these remaining elements are taken in 
account. Let us define by Ca (respectively Cx) the num- 
ber of blocks containing only elements of species A {X 
respectively) which will be created with the (respec- 
tively rx) remaining elements. We propose the following 
alternatives for the choice of Ca and Cx ■ 

1. liftyi Given rA and rx, we build as many blocks as 
possible, i.e., 



C 



X 



rx 



1 



and Ca 



rA 



1 



Then, most of the Ca blocks of species A will have 
length La + 1; others will be longer because of the 
remainder of the division between rA and La + 1- 
The same is true for the clusters of species X. Once 
these blocks have been created we have only to de- 
cide how to place them on the lattice. Having most 
of the blocks of equal length Lx + 1 may introduce 
artifacts on the distribution of clusters resulting 
from the lifting. To partially increase the random- 
ness of such distribution we modify the equations 
for Cx and C a by taking them to a random integer 
number picked by a uniform distribution between 



rx 



rx 



L 



X 



2' L 



X 



1 



and 



rA 



rA 



La + 2 ' La + 1 



(20) 



In such way the length of the blocks is distributed 
more randomly. We show in the next section how 
this choice for Cx and Ca still introduces artificial 
oscillations in the distribution of clusters. 

2. lift B Instead of building as many blocks as possible, 
we might generate a more realistic number of them. 
To this end, we derive an upper and a lower limit for 
Ca and Cx , and we additionally use the constraint 
that the total number of clusters for both species 
is equal, i.e.. 



CA + mA=^Cx + mx- 



(21) 



Then, the maximum number of clusters for species 
X can be derived to be 



Cm ax 
X 



-\- max(m^ — mx, 0) — max(TOx — mA, 0)^ , 



( 


rx 




rA 




[Lx + 1\ 




La + 1_ 



whereas for the minimum we have 



C™" = max((5r^ , tua - mx), 



(22) 



where Sr^ = 1 if ?'x > and otherwise; Srx 
takes into account the unlikely event where all the 
elements X are contained into Bx- Similar formu- 
las can be derived for C™^'^ and C^"". We then 



pick a number of clusters Cx randomly, based on 
a uniform distribution between the corresponding 
minimum and maximum, and infer the number of 



clusters Ca from (21 1 



3. liftc In lifts, the lower limits C™'^ can be very 
conservative, in fact this limit refers to a configura- 
tion where all the elements of rx (respectively r^) 
are packed together in a single long cluster, such 
that the resulting lifting operator artificially cre- 
ates a low number of clusters. This corresponds 
to clusters which are artificially long. We there- 
fore also investigate a lifting in which the minimum 
number of clusters is based on the derivation of a 
"reasonable" maximal cluster length, based on the 
specific values of rx and r^, from which the min- 
imal number of clusters C™" then is derived, the 
argument based on the probabilistic reasoning. In 
the next section, we write an explicit estimate for 
the example |I.1[ 

Once the numbers of clusters have been determined, 
and the corresponding individual clusters have been cre- 
ated, they need to be positioned on the lattice. Let us 
denote by (respectively B^) the sets containing the 
clusters built out of the remainders of each species. Re- 
call that, given the assumption that mx > rriA, we still 
need to place, besides the clusters corresponding to the 
remainders, mx — mA clusters of species X correspond- 
ing to the macroscopic state variables. To this end, we 
continue as before, only now we alternatingly pick an ele- 
ment from Bx and i?^. After that step, for lift a the re- 
maining := CA + Cx^{mx~rnA) clusters in B^UB^ 
are randomly picked and appended to the configuration; 
for lifts and liftc instead we have an equal number of 
elements remaining in B^ and B^ so it is sufficient to 
pick in the alternate way from the two sets. 



Note that from the cluster definition (12 1 when recon- 
structing we impose adjacent clusters to superimpose at 
the head and the tail. The case where mx > mA can be 
handled identically. 

Remark III. 1 (Lifting based on coverage) . When choos- 
ing La = Lx = 0, the only macroscopic state variable is 
the coverage a, i.e. M — {a}. In that case, the lifting 
generates a configuration cr in which the position of the 
N X a elements of species A is completely random. This 
method is labeled as liftfl. 



IV. NUMERICAL EXPERIMENTS 

In this section, we illustrate the use and the efhciency 
of the different lifting operators defined above. We ap- 
ply them on the two examples introduced in Section [ij 
analyze and quantify the deviations from the original mi- 
croscopic simulations and discuss the reasons behind the 
observed discrepancies. 
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FIG. 1. Results for the microscopic dynamics for K — K\ 
and K — K2 reaction rates. 



A. Example [TT] 



We study the system ([2| with two different reaction 
rates sets: Ki — [ki = 2,fc2 ~ l,fc3 = 0.1, /c4 — 0.01] 
and K2 = [ki = l,k2 = 2,/c3 = 0.01,^4 = 0.1]. In both 
cases, the single-element reactions are slower than the 
reactions including nearest-neighbor interactions. The 
choice of these specific parameter sets is motivated by 
the differences in the resulting equilibrium values. For 
the reaction rates Ki, the concentration at equilibrium 
is Oeq — 0.9; for the reaction rates K2, is acq — 0.35, 
see Figure [l] Unless otherwise stated, all simulations 
are performed on a lattice with N = 2000 sites for a 
time interval T = [0,400]. We show in the next sections 
how the different equilibrium values effect the quality of 
the various liftings. We first assess the accuracy of the 
lifting procedure lift^j whose macroscopic description is 
solely based on coverages (Section |IV A 1 1. We then in- 
troduce the augmented macroscopic variables set (Sec- 
tion IV A 2 ) and perform a numerical test using each of 



the discussed lifting strategies (Sections IV A 3 TV A 5 1. 
The test consists of simulating using the coarse time- 
stepper using the minimal possible time step. Hence, 
after every microscopic time step, we first restrict to the 
selected macroscopic state variables, after which we lift 
back to a microscopic configuration. As a result, the 
observed macroscopic time evolution will mimic the nu- 
merical closure approximation induced by the lifting op- 
erator. Finally, we analyze the dynamics generated by 
these algorithms for larger and larger coarse time-steps 
(i.e., increasing the time between liftings), we will discuss 
the effects and limitations of the "healing principle" in 



this context (Section IV A 6). 



1. Coverage is insufFicient as macroscopic state variable 

We first present results using the coverage a as the only 
macroscopic state variable, from which we lift using the 
procedure lift^, see Remark [lILl[ The numerical results, 
shown in Figure [2] indicate a significant deviation of the 
dynamics, as well as of the obtained equilibrium. In fact, 
in this case the lifting simply leads to the corresponding 
mean field dynamics. This result is not surprising per 
se: when taking into account only the coverages, and 
randomly placing the elements of species A, one is in 
fact stirring the system, so spatial correlations are lost, 
and the system effectively satisfies the requirements to 
display mean field behavior. 

In the literatur^^ it has been argued that errors in 
the macroscopic evolution that are induced by an incor- 
rect microscopic initial condition may be able to "heal" 
by taking a larger coarse time step St. Therefore, we 
enlarge the size of the coarse time step (to allow for heal- 
ing). Figure[2]shows that this leads to an improved agree- 
ment. However, for our system, the healing effect does 
not suffice for good macroscopic evolution. In our numer- 
ical tests we see that performing few hundreds of lifting 
procedures in a long (hundreds of thousands steps) simu- 
lation is sufficient to spoil the microscopic dynamics and 
to modify the equilibrium values. 

Obviously, using only the coverages as macroscopic 
variables does not seem to be, generally speaking, a good 
starting point for lattice reactive systems. 



2. Selection of macroscopic state variables 

Given the results from the previous section, we aug- 
ment the set of macroscopic state variables which in turn 
induces the necessity to choose a more sophisticated lift- 
ing procedure (as discussed above). As a first sanity 
check for the methodology, we initially perform a sim- 
ulation with the coarse time-stepper using a very high 
number of macroscopic state variables, i.e., Lx — La — 
N/20. Obviously, in this case we perform almost no re- 
duction, in the sense that only a few microscopic config- 
urations are consistent with the considered macroscopic 
state. Figure [3] shows the results for both selected sets 
of reaction rates in the case of for lifts, equivalent re- 
sults were obtained with the other lifting procedures but 
are not shown here. As expected, the agreement with the 
full microscopic simulations is almost perfect in this case. 
Our objective here is however to find the simplest (and 
most effective) algorithm that correctly reproduces these 
simulations. In this context, one will want to determine 
and use the smallest possible set of macroscopic variables 
for the corresponding coarse time-stepper. We propose 
to base this choice on the previously introduced dynam- 
ical equivalence, which in our sense is a good measure of 
the instaneous "reactivity" of a system. 

Remember that we have considered two microscopic 
configurations to be equivalent if in both configurations. 
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FIG. 2. Failure of healing, results for lift_R for examples |I.1| 
with K — Ki (up) and with K = K2 (down). The results 
refer to runs with 1000 microscopic time steps as coarse time 
stepper between every lifting/restriction procedure (Hi) and 
with 10000 microscopic time steps {H2)- Mean Field refers to 
Q. The total simulation time was 300000 steps. 



every reaction in the system is equally likely, i.e., in both 
configurations, there is an equal amount of reactive ele- 
mentary blocks, see equation (16). For the chemical re- 
actions ([2]), we observe that the first reaction introduces 
a dependence on the neighboring sites, while the second 
one does not. Consequently, when choosing La and Lx 
we focus on preserving the number of possible reactants 
of the first reaction, as the possible reactants for the sec- 
ond reaction are already determined by the coverages. 

Let us first decide upon La- In particular, we note that 
isolated As are the only possible elements which may be 
modified by the first reaction. Therefore, we choose La — 
1. The choice of Lx is more difficult, as one can never 
preserve the number of elements X which are surrounded 
by Xs, unless Lx is very large. However, we anticipate 
(and will show later on) that by keeping La = 1, while 
increasing Lx, we should be able to get to the smallest 



FIG. 3. Comparison of microscopic (continuous line) and 
restriction/lifting (dashed line) runs with La ~ 100 = Lx = 
100. The results shown are for reaction rate Ki (up) and K2 
(down). The runs refer to lifts, but equivalent results are 
obtained for the other lifting methods. 



set of macroscopic variables for which the coarse time- 
stepper recovers the macroscopic dynamics of the original 
microscopic system. For small values of Lx , we will show 
that the quality of the recovered macroscopic dynamics 
depends on the specifics of the lifting operator. 



3. Lifting technique liftA 

We now proceed by performing a simulation using the 
coarse time-stepper with the lifting lift^, using La — 1 
and different values of Lx- The results for both sets of 
reaction rates {Ki and K2) are shown in Figure |4] for the 
smallest possible coarse time step. For the set of reac- 
tion rates Ki , when increasing Lx , we quickly obtain an 
excellent agreement between the numerical closure ap- 
proximation and the fully microscopic simulation. For 
the set of reaction rates K2, however, we see that for 
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FIG. 4. Restriction/lifting results for liftA for both sets of 
reaction rates K = (up) and K = K2 (down) . Comparison 
between different choices of Lx ■ 



lift^, starting from the macroscopic state variables M\ 
and (1 < ^ < Lx), for different values of Lx, and 
average again the number of clusters of length 1 < ^ < 
N over the ensemble of microscopic configurations (see 
Figure [5]). 

For the original microdynamics, this distribution de- 
cays exponentially with the length of the clusters. Such 
rapid decrease can be explained by observing the reaction 
laws ([2|. Any long sequence of X is likely to be broken 
in two smaller ones by both reactions, while sequences of 
As can be interrupted only by the second reaction, so the 
dynamics greatly favors short clusters of Xs rather than 
longer ones. This feature is not preserved by the lifting 
procedure; as can be seen in Figure [5j for the lifted dis- 
tributions, longer clusters are much more likely than for 
the micro-dynamics. Moreover, we also observe oscilla- 
tions in the different cluster populations depending on 
the choice of Lx which determine the size of the blocks 
that are contained in the remainder B^^. 

While appealing due to its simplicity, the lift^ proce- 
dure thus suffers from some major drawbacks that can 
be related to its inability to correctly preserve the dis- 
tribution of clusters size. It is relevant at this stage to 
identify the reasons behind this inaccuracy. Recall that 
in liftyi, the idea is to maximize the number of clusters 
of species A and X created with the remainders, but 
that no special attention is being paid to the relative val- 
ues of Ca and Cx- As a result, in the final step of the 
lifting, multiple building blocks of the same species can 
be selected, resulting in longer clusters. To analyze the 
lifted distributions, we need to estimate the probability 
of picking in sequence k blocks corresponding to species 
X from U B\ (see end of Section III ) . For small val- 



ues of fc, relative to the total number of available blocks, 
i.e., k vT , this probability can be approximated by a 
simple multiplication of the probability of picking first a 
block of A then k blocks X and then another block X: 



similar values of Lx, the macroscopic evolution behaves 
completely differently (both dynamically and in equilib- 
rium), and we need to increase Lx quite substantially 
before getting an agreement. 

This difference can be rationalized by analyzing the 
distribution of cluster size before and after the lifting. 
To this end, we start from an equilibrium microscopic 
configuration for a ring of length N = 2000 and count 
the number of clusters and {I = 1, • ■ • , N) via 



equation ( 12 ) at each time step during a microscopic sim- 
ulation, for a total time of T = 100 and average over 
time. 



cq. 



cq. 



where the subscript eq emphasizes the fact that the time 
average is taken over a simulation that is in equilibrium 
for all values of t. The resulting (normalized) histogram 
is plotted in Figure [s] (red starred line). Next, to obtain 
a probability distribution for the coarse time stepper, we 
repeat the lifting procedure 10000 times using method 



where px 



Vk^p'xil'Pxfn^, (23) 
^ is the probability of getting a block 



of Xs out of the whole set U B^- Using (|23|), and 



considering that the cluster population is simply kVk, 
we see that if px — 0.5, i.e. we have approximately 
equally many building blocks for both species, then we 
get almost a plateau in the populations for clusters of 
length slightly larger than Lx ■ In Figure [5] we see such 
a plateau when Lx = 0. Remember that in this case, 
all elements of species X are contained in the remainder 
B^, as the set Bx is empty. When choosing Lx > 0, the 
behavior is different. In that case, we observe oscillations, 
which can be explained by noting that clusters are built 
by concatenation of building blocks, and hence, multiples 
of Lx + 1 are more likely to occur. 

Increasing Lx gives better results, because the number 
of building blocks of species X in the remainder dimin- 
ishes, and, moreover, the size of the those blocks also 
increases. 



11 



10 ■ 



10 ' 



10" 



10' 



10" 



10" 



micro 
•-• liftA L_X=0 
liftA L_X=3 
liftA L X=5 




10" 
10^ 
10^ 
10^ 
10' 
10' 

lo' 



i 6 8 10 12 

cluster size 

X cluster population 



14 



micro 
liftA L_X=0 
liftA L_X=3 
♦ liftA L X=5 




10 15 20 
cluster size 



25 



30 



FIG. 5. Comparison of cluster populations for species X for 
liftA at equilibrium in a semilogarithmic scale. The red stars 
correspond to the average over a microscopic run at equilib- 
rium for K — Ki (up), and K = K2 (down). The other data 
correspond to the average of 10000 liftings, given the values 
of the macroscopic set at equilibrium. 



4. Lifting technique lifts 

We now repeat the previous experiment for the lifting 
operator lifts, thus starting with the minimum coarse 
time step. For the reaction rates Ki, we observe a much 
better agreement with the microscopic data for small val- 
ues of Lx- Choosing Lx = 1 already gives excellent 
agreement, see Figure |6] (left). For K2, instead, the re- 
sults are at first glance somewhat surprising. In fact, 
we get a very good agreement between the coarse time- 
stepper and the microscopic simulation by using only 
Lx = and La = 1- However, when we add extra 
macroscopic variables, the results initially worsen. 

The main reason for this unanticipated behavior may 
be found by analyzing the dependence of the minimal and 
maximal possible number of clusters in the remainders 
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FIG. 6. Restriction/lifting results for lifts for both sets of 
reaction rates K\ (up) and K2 (down); the figure on the right 
is zoomed to the equilibrium region. 



and La- We observe that a number tua — 'nix of blocks 
of species A may be left over after placing all elements 
contained in Bx- These elements need to be placed in 
the final step of the lifting, together with the remain- 
ders of species X. For fixed L^, the number niA — mx 
reaches its maximum at Lx = 0, so C^'" is maximal 
for Lx = 0. Simultaneously, because of the fact that 
_B^. then contains all the elements of species X present 
in the system, Cx'^^ is maximal for Lx = 0. A high 
value for (7™'"'™'^'^ implies that the elements of species 
X will be distributed in a large number of short clusters, 
which corresponds (by coincidence) to the structure of 
the microscopic realizations. 

When Lx > 0, the quantity tua — mx decreases, and 
so does C^'", until it reaches its minimum value of 1, 
corresponding to a case in which all remaining elements 
of species X are contained in a single long cluster. (Note 
that, theoretically, one may also encounter situations in 
which the set is empty, corresponding to Cx =0.) 
It is then that we observe the largest error for the lifting. 
After this peak, the lifting performs better, as the num- 
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FIG. 7. Effect of increasing Lx using lifts and liftc when 
K — K2- The error represent the time average of the square 
difference between the a coverage data coming from a micro- 
scopic run and a lifting/restriction one with initial conditions 
at equilibrium. 



ber of elements in B^^ decreases, leaving less room for 
artifacts in the evolution (Figure |6] and Figur^ ) . As for 
lift^, we have evaluated the X cluster population at the 
equilibrium for K = K2 and the results, shown in Fig- 
ure [s] (left), are in line with the above observations. We 
note that the decay in cluster populations as a function 
of cluster length now becomes non-oscillatory. Neverthe- 
less, this decay is still slow compared to that observed in 
the full microscopic simulation. 

To summarize, when increasing the number of macro- 
scopic variables Lx , we have a competition between two 
opposite effects. On the one hand, as expected, tak- 
ing into account more variables reduces the degrees of 
freedom that are affected by the lifting, increasing the 
accuracy. On the other hand, we observe that, for in- 
termediate values of Lx, a bias is introduced (via C^'" 
in this case) in the distribution of the reconstructed el- 
ements. Depending on the relative importance of these 
effects, convergence to the microscopic dynamics may be 
non-monotonic. Note that this does not explain the very 
good agreement for the case Lx = 0, which is coinci- 
dental, and appears to be caused by the fact that lifting 
errors that artificially generate clusters that are too short 
do not affect the macroscopic dynamics. 



5. Lifting technique liftc 

Finally, we repeat the experiment for the lifting oper- 
ator liftc. The idea is here to use an estimate of the 
cluster size distribution to find an appropriate estimate 
for the minimum number of clusters: we may for example 
use ( 23 ) to get an estimate for the maximum length, from 

k 



10" 
10^ 
10^ 



-5 10 



X cluster population 



10 



10 " 



10 



10" 



10" 



lO' 



10" 



10"' 



10" 



10"' 



micro 
liftB L_X=0 
liftB L_X=5 
liftB L_X=8 
lifts L X = 10 




10 



15 

cluster size 



20 



25 



30 




micro 
liftc L_X=5 
liftc L_X=8 
liftc L X=10 



10 
cluster size 



15 



20 



FIG. 8. Comparison of X cluster population for lifts (up) 
and liftc (down) at equilibrium in a semilogarithmic scale, 
the micro data refer to the average over a microscopic run at 
equilibrium for K2 , while the lifting are the repetition of 10000 
liftings given the values of the macroscopic set at equilibrium. 
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FIG. 9. Restriction/lifting results for liftc and K — K2 
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with k such that 



pi(i-Px)V 



< e, 



(24) 



for a chosen small value of e, explicitly excluding config- 
urations that are very unlikely from a probabilistic point 
of view. (Here, we chose e = 0.03.) In our numerical 
experiments, the macroscopic simulation results for the 
reaction rates Ki were very similar to those obtained 
with lifts . This is not surprising, as the lifting tech- 
nique is similar in the sense that it also picks the number 
of clusters homegeneously between two imposed values, 
the only difference between these two schemes being in 
the estimate of the minimal number of clusters. We nev- 
ertheless notice that, for K2, the results improve slightly. 
Considering again the distribution of cluster sizes (see 
Figure [S] (right)), we observe that some oscillations on 
the cluster populations are still present, but overall the 
distributions obtained with the liftings are closer to the 
data coming from the microscopic dynamics thanks to a 
better estimate for one of the boundaries in the number 
of clusters. 



6. Effects of healing 

As previously mentioned we have also checked the ef- 
fect of the healing on the liftings lift g and liftc for those 
macroscopic variables sets which present the largest de- 
viation from the microscopic dynamics, see Figure |10[ 
Healing is expected to occur whenever a sufficiently long 
coarse time step is chosen. The two lifting procedures re- 
act in a quite different way to the increase of the coarse 
time step. We observe that even considering very large 
healing time, with lifts we still do not recover the proper 
dynamics, instead for liftp we observe a very good agree- 
ment using a healing time of the order of 1/300 of the to- 
tal time. The reason for such result may again be founded 
on the fact that with liftc the error on the cluster dis- 
tribution for species X are smaller than with lifts and 
such small error are corrected by short healing. 



B. Example [O] 

Let us now turn to the simpler Example |I.2[ in which 
only the nonlinear reaction is present. A direct micro- 
scopic simulation shows that a 1-dimensional ring, ini- 
tially completely covered by X evolves towards a macro- 
scopic steady state with Ocq — 0.28. If, instead, we look 
at the solutions of the approximating mean field equa- 
tions ( 11 ), we obtain a steady state with equal concentra- 

— — 0.5. 



tions of both species, and consequently a^q ~= x, 
Obviously, a more sophisticated approach is here also 
necessary. 

We first consider the lifting liitfj, see Remark |III.1[ 
which only takes into account the coverage of species A, 
and perform the same numerical experiment as for Ex- 
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FIG. 10. Restriction/lifting results for lifts (up) and liftc 
(down) for K — K2 with healing. In both the cases we com- 
pare the lifting procedure with and without healing (of 1000 
microscopic time steps) The total simulation time is T = 400, 
1000 steps correspond roughly to a time of 1. For the param- 
eters we have A'^ = 3000 and Lx ~ 5. 



ample 
urelTT 



[LTJ see Section pV A| The results are shown in Fig- 
We observe that, as expected, liftn has a = 0.5 as 
steady state, demonstrating that this lifting induces the 
same closure approximation as the mean field equations. 

A significant improvement can be obtained by making 
the following observation: if we start from a uniform X 
coverage, then any element of species A appearing in the 
ring is always necessarily surrounded by two elements of 
species X. Therefore, any cluster longer than a isolated 
A is physically impossible. As a consequence, when we 
calculate the coverage, we are already calculating Al\, 
and we can easily take {a, M\} as the macroscopic vari- 
ables set, instead of the sole coverage a. Note that, for 
this choice of macroscopic state variables, all lifting pro- 
cedures introduced in Section |III| become very similar. 
Indeed, for any value of Lx, the remainder of elements 
of species A that are not accounted for by the macro- 
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FIG. 11. Comparison of micro dynamics with lift_R and 
lift A, both using La = 0, Lx ~ 0. 



scopic state variables is empty (_B^ = 0). The only choice 
that needs to be made is the division of the elements of 
into clusters of different length that are necessary to 
separate all elements of species A. Here, we choose to 
simply separate each two elements of species A by a sin- 
gle cluster of species X of length Lx + 1- The remaining 
elements of species X are then inserted randomly next 
to an element of species A. We call the resulting lifting 
operator lift a • 

The macroscopic dynamics that results from the coarse 
time-stepper is also depicted in Figure [TT] The result 
illustrates that it is enough to take out the physically 
impossible states to get an excellent agreement between 
the microscopic and the coarse-grained dynamics. We 
reiterate that both liftings essentially start from the same 
macroscopic information, namely the number of elements 
of species A present in the system; the only difference is 
that, for liftyi, we have used some additional information 
about the specific features of the microscopic dynamics, 
namely the fact that a = M\. 



V. DISCUSSION AND OUTLOOK 

In this paper, we have presented and studied an ap- 
proach to implement numerical closure strategies, aimed 
at reproducing the dynamical behavior of stochastic re- 
active systems taking place on low-dimensional supports 
(here, on a lattice). Such closure is desired in the frame- 
work of multi-scale simulations, where one needs to ac- 
celerate the original, fully microscopic simulations. It 
is also particularly relevant for performing "upper-level" 
tasks on the system, such as bifurcation analysis or con- 
trol design. Performing these tasks is often challenging 
in the case of low-dimensional systems, as it is not possi- 
ble to use macroscopic approaches like the mean field 
approximation to have a reliable macroscopic descrip- 



tion of the dynamics. The multi-scale strategies we used 
here are based on the idea of constructing a coarse time- 
stepper for a suitable, preferentially small set of macro- 
scopic state variables. Microscopic simulations are then 
only called for to estimate the (local) time derivative of 
such macroscopic variables. The main difficulty in such 
an approach lies in the definition of a suitable lifting oper- 
ator that maps the macroscopic state variables to a given 
microscopic configuration. The investigation of different 
choices for the lifting was exactly the focus of the present 
paper. In particular, the main goal was to investigate 
which lifting operator recovers the macroscopic dynam- 
ics of the system with a minimal number of macroscopic 
variables. 

We thus investigated several approaches using the 
number of clusters of particles as macroscopic vari- 
ables for the system. This choice allows for a rela- 
tively straightforward implementation of a lifting oper- 
ator, which then consists in first placing such elementary 
blocks on the lattice, and then filling the rest of the sys- 
tem with larger clusters constructed with the remaining 
particles. We proposed and analyzed several lifting pro- 
cedures that differ in the way they estimate the size and 
the amount of such larger blocks. We performed several 
numerical experiments with these different approaches, 
which demonstrated how few variables might be sufficient 
to obtain a good agreement between the macroscopic and 
the microscopic dynamics. In particular, with the same 
set of macroscopic variables we may obtain significantly 
different results using different lifting strategies. 

We specifically identified two main possible ways to 
perform the lifting. There exist "simple" strategies where 
it is assumed that nothing is known on the microscopic 
dynamics and the population of larger clusters is esti- 
mated in a "rough" way. There are also "complex" 
strategies, where more information on the microscopic 
behavior is taken into accountn such as a simulation- 
based estimate of the cluster size distribution. For the 
"simple" lifting we note that the accuracy always in- 
creases when increasing the size of the macroscopic vari- 
ables set, but that at the same time a high number of 
them is needed to get a satisfactory accuracy. In the 
case of "complex" liftings, we observe a general better 
agreement with the microscopic data points, but it is 
possible that convergence to the microscopic dynamics 
is non-monotonic as a function of the number of macro- 
scopic state variables. Within these possible approaches 
we have also studied the effect of applying a long coarse 
time step, allowing the system to "heal" between different 
lifting procedures. We showed how, for some liftings, re- 
lying on healing is not sufficient to correct errors coming 
from the reconstruction, but also that in some cases the 
healing substantially improve the macroscopic results. 

Although we could identify efficient lifting strategies 
for the case hereby considered, it is not entirely clear how 
to best choose such strategy in the most general case. It 
seems that there exists a trade-off between the number of 
macroscopic state variables that one retains, versus the 
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amount of effort that is put in the initiahzation of the 
remaining degrees of freedom of the system. An inter- 
esting possibihty would be to eonstruet an approximate 
distribution of large clusters, based on the most pop- 
ular analytical closure schemes (such as the Kirkwood 
approximation, the Ursell expansion, etc.) and then rely 
on the healing of the system during the coarse time step 
to correct the (hopefully) small errors introduced in this 
way. This approach would in addition allow one to for- 
malize better the action of the lifting strategies and help 
choosing the best procedure, based on the known physical 
properties of the system. The efficiency of a particular 
choice of numerical closure is not only relevant for simu- 
lation purposes, but also give important information on 
what closure relations make sense in the system under 
consideration. It is our opinion that the establishment 
of such a connection would open the way to a system- 
atic numerical investigation of the level of accuracy of 
existing possible closures in stochastic systems. 
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